function D = chebop(N, xmin, xmax)
  x = chebpoints(N, xmin, xmax)';
  factor = 2./(xmax-xmin);
  N = N - 1;
  c = [2; ones(N-1, 1); 2].*(-1).^(0:N)';
  X = repmat(x, 1, N+1);
  dX = X - X';
  D = (c*(1./c)')./(dX+(eye(N+1)));
  D = D - diag(sum(D'));
  D = D;
end

